// includes, system
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <assert.h>
#include <float.h>

// includes, project
#include <cutil.h>

// includes, kernels
#include <pdist_kernel.cu>
#include "pdist.h"
#include "linkage.h"
#include "cluster.h"

/**
 * Initialize CUDA.
 *
 * Returns:
 * 0 - CUDA initialized properly
 * 1 - No device found
 * 2 - Device does not support CUDA
 */
extern "C" __declspec(dllexport)
int init_cuda()
{
	// taken from CUT_DEVICE_INIT
	int deviceCount;
	int dev;

	CUDA_SAFE_CALL_NO_SYNC(cudaGetDeviceCount(&deviceCount));
	if (deviceCount == 0) {
		return 1;
	}
	for (dev = 0; dev < deviceCount; ++dev) {
		cudaDeviceProp deviceProp;
		CUDA_SAFE_CALL_NO_SYNC(cudaGetDeviceProperties(&deviceProp, dev));
		if (deviceProp.major >= 1)
			break;
	}
	if (dev == deviceCount) {
		return 2;
	}
	else {
		CUDA_SAFE_CALL(cudaSetDevice(dev));
	}

	return 0;
}

/**
 * Returns the lower triangle of a square matrix.
 *
 * output can be the same address as input.
 */
extern "C" __declspec(dllexport)
void lower_triangle(float *output, const float *input, unsigned int size_y)
{
	unsigned int x, y;
	unsigned int counter = 0;

	for (y = 0; y < size_y; y++)
	{
		for (x = y+1; x < size_y; x++)
		{
			output[counter] = input[y*size_y + x];
			counter++;
		}
	}
}

/**
 * Compute the pairwise Euclidian distance among rows of the input matrix.
 *
 * Be sure to run init_cuda() before calling this function.
 */
extern "C" __declspec(dllexport)
void run_pdist(float *output, float *input, unsigned int size_x, unsigned int size_y)
{
	const unsigned int input_size = sizeof(float) * size_x * size_y;
	const unsigned int output_size = sizeof(float) * size_y * size_y;
	float *gpu_input, *gpu_output;

	CUDA_SAFE_CALL( cudaMalloc( (void**) &gpu_input, input_size) );
	CUDA_SAFE_CALL( cudaMalloc( (void**) &gpu_output, output_size) );
	CUDA_SAFE_CALL( cudaMemcpy(gpu_input, input, input_size, cudaMemcpyHostToDevice) );

	dim3 grid(1, (size_y%BLOCK_DIM == 0 ? 0 : 1) + size_y / BLOCK_DIM, 1);
	dim3 threads(1, BLOCK_DIM, 1);
	pdist_gpu<<< grid, threads, sizeof(float)*size_x >>>(gpu_output, gpu_input, size_x, size_y);
	cudaThreadSynchronize();
	CUDA_SAFE_CALL( cudaMemcpy(output, gpu_output, output_size, cudaMemcpyDeviceToHost) );
	CUT_CHECK_ERROR("pdist: Kernel execution failed");

	CUDA_SAFE_CALL( cudaFree(gpu_input) );
	CUDA_SAFE_CALL( cudaFree(gpu_output) );
}

/**
 * Parse a file into a matrix.
 * 
 * This function does not deal well with errors. It should only be used for
 * testing purposes where the input is guaranteed to be in the correct format.
 */
static void parse_file(const char *filename, float *array, int size_x, int size_y, int transpose)
{
	FILE *f;
	char *line, *orig;
	int x = 0, y = 0;

	line = (char *)malloc(sizeof(char) * 1024);
	orig = line;

	f = fopen(filename, "r");
	//assert(f);
	while (fgets(orig, 1024, f) && y < size_y)
	{
		line = orig;
		x = 0;

		while (line)
		{
			char *old = line;
			line = strstr(line, " ");
			if (line)
			{
				*line = '\0';
				line++;
			}

			if (transpose)
				array[y + x*size_y] = (float)atof(old);
			else
				array[y*size_x + x] = (float)atof(old);

			x++;
		}

		y++;
	}

	fclose(f);
	free(orig);
}

struct test {
	char file[128];
	int x;
	int y;
};

/**
 * Run tests to determine validity of the calculations.
 */
void runAllTests(int argc, char **argv)
{
	const int num_tests = 8;

	const struct test tests[num_tests] = {
		{"y10x10", 10, 10},
		{"y100x20", 20, 100},
		{"y1000x20", 20, 1000},
		{"y3000x16", 16, 3000},
		{"y4000x8", 8, 4000},
		{"y5000x7", 7, 5000},
		{"y6000x32", 32, 6000},
		{"y7680x16", 16, 7680},
	};

	char filename[128];
	float *output, *input, *reference, *result;
	int i;

	for (i = 0; i < num_tests; i++)
	{
		int size_x = tests[i].x, size_y = tests[i].y;
		printf("Testing %s... ", tests[i].file);

		input = (float*)malloc(sizeof(float)*size_x*size_y);
		memset(input, 255, sizeof(float)*size_x*size_y);
		output = (float*)malloc(sizeof(float)*size_y*size_y);
		memset(output, 0, sizeof(float)*size_y*size_y);
		reference = (float*)malloc(sizeof(float)*3*(size_y));
		memset(reference, 0, sizeof(float)*3*size_y);
		result = (float*)malloc(sizeof(float)*3*(size_y));
		memset(result, 0, sizeof(float)*3*size_y);

		sprintf(filename, "tests/%s.in", tests[i].file);
		parse_file(filename, input, size_x, size_y, 1);
		sprintf(filename, "tests/%s.out", tests[i].file);
		parse_file(filename, reference, 3, size_y-1, 0);

		run_pdist(output, input, size_x, size_y);
		lower_triangle(output, output, size_y);
		run_linkage(result, output, size_y);

		if (compare_matrix_linkage(result, reference, 3, size_y-1) == 0)
			printf("OK\n");
		else
			printf("FAILED!\n");

		free(input);
		free(output);
		free(reference);
		free(result);
	}
}

int main( int argc, char** argv) 
{
	CUT_DEVICE_INIT();

	runAllTests(argc, argv);

	return(EXIT_SUCCESS);
}

